#!/usr/bin/env python
 
'''
pso.py
 
A simple implementation of the Particle Swarm Optimisation Algorithm.
Uses Numpy for matrix operations. 
Pradeep Gowda 2009-03-16

modified by 
Valerio Amodeo and Marc Baro
2013-11-12

'''
 
from numpy import array
import scipy.spatial.distance as Ssd
from random import random, randint, uniform
from math import sqrt, exp, fabs
from sys import stdout
import Gnuplot
 
#constants and functions parameters 
max_iterations = 1000
goalFitness = 0.001
'''
#population_size = 100
c1 = random()*2
c2 = random()*2
globalBest=None
particles = []
fitnessVector = []
i = 0
'''

class Particle:
    pass
        
class Pso:

	def __init__(self, pop, dim, c1, c2):
		self.i=0
		self.packages = []
		self.population_size = pop
		self.dimensions = dim 
		self.fitnessVector = []
		self.particles =[]
		self.valueWeight = None
		self.c1 = c1
		self.c2 = c2 
		self.globalBest=None  

	@staticmethod
	def circleFitness(param):
		fun = 0
		para = param*10 #removed this to test
		#f(u1,..,ux) = ((u1*u1)+ ...  +(ux*ux) 
		for i in para: fun += i*i
		return fun

	'''def knapsackFitness(self, param):
		value = 0
		weight = 0
		for i in range(len(param)):
			#print 'i = ', i, '  pack = ', self.packages[i] 
			weight += (0 if param[i] < 0.5 else 1)*self.packages[i][1] #weight			
			value += (0 if param[i] < 0.5 else 1)*self.packages[i][0] # value
		self.valueWeight = (value, weight)
		if weight > 1000:
			return -1
		else:
			return 100*(1/(1+fabs(1000-weight))) + 50*value
	'''
	def knapsackFitness(self, param):
		value = 0
		weight = 0
		for i in range(len(param)):
			#print 'i = ', i, '  pack = ', self.packages[i] 
			weight += (0 if param[i] < 0.5 else 1)*self.packages[i][1] #weight			
			value += (0 if param[i] < 0.5 else 1)*self.packages[i][0] # value
		self.valueWeight = (value, weight)
		fun = value - 0.1*fabs(0 if 0 < (1000-weight) else (1000-weight))
		#print fun
		return fun
	
	def initParticles(self, maxValue):		
		for i in range(self.population_size):
			p = Particle()
			if maxValue:
				p.params = array([random() for i in range(dimensions)])
				p.fitness = 999999999.0
			else:
				p.params = array([(random() if random()<0.1 else 0) for i in range(dimensions)])
				p.fitness = 0
			p.v = 0.0
			p.best = p.params
			self.particles.append(p)
			
	def volumeKnapsackFitness(self, param):
		fun = 0
		#f(u1,..,ux) = u1*p1.val + u2*p2.val + ... + ux*px.val
		weight = 0
		volume = 0
		for i in param: 
			weight += (0 if i < 0.5 else 1)*packages[i][1] #weight
			volume += (0 if i < 0.5 else 1)*packages[i][2] #volume
			if weight > 1000 or volume > 250:
				fun -= (0 if i < 0.5 else 1)*packages[i][0] #value
			elif weight <= 1000 and volume <= 250:
				fun += (0 if i < 0.5 else 1)*packages[i][0] #value
		return fun
	
	def initPackages(self):
		f = open('pso-packages')
		for line in f:
			string = line.split(',')
			self.packages += [(float(string[0]), float(string[1]))]
	
	def initPackagesVolume(self):
		f = open('pso-packages')
		for line in f:
			string = line.split(',')
			self.packages += [(float(string[0]), float(string[1])), uniform(1.0,5.0)]
	
	@staticmethod
	def weight(t, lam=0.00091):
		return exp(-lam * t)
	
	def findBestNeighbour(self, particle):
		neighbours = []
		minDistance = 0
		#select a random particle
		current = self.particles[randint(0, (len(self.particles)-1))]
		#calculate the euclidean distance
		minDistance = Ssd.euclidean(particle.params, current.params)
		while len(neighbours)<3:
			for p in self.particles:
				if p not in neighbours:
					distance = Ssd.euclidean(particle.params, p.params)
					if distance < minDistance:
						current = p
						minDistance = distance
			neighbours += [current]
			while True:
				current = self.particles[randint(0, (len(self.particles)-1))]
				if current not in neighbours:
					break
			distance = Ssd.euclidean(particle.params, current.params)
		current  = neighbours[0]
		for p in neighbours:
			if p.fitness < current.fitness:
				current = p
		return current
		
	def particleSwarm1(self):
		v = 0
		self.globalBest = self.particles[randint(0, (len(self.particles)-1))]
		while self.i < max_iterations :
			for p in self.particles:
				fitness = Pso.circleFitness(p.params)
				if fitness < p.fitness:
					p.fitness = fitness
					p.best = p.params
				if fitness < self.globalBest.fitness:
					self.globalBest = p
				# v(t+1) = v(t) + (c1 * r1 * (p(t) - x(t)) + (c2 * r2 * (g(t) - x(t))) 
				v = p.v + c1 * random() * (p.best - p.params) \
						+ c2 * random() * (self.globalBest.params - p.params)
				# x(t+1) = x(t) + v(t+1)
				p.params = p.params + v
			#saving the best Global Fitness per cicle for the plot
			self.fitnessVector += [(self.i+1, self.globalBest.fitness)]
			self.i  += 1
			stdout.write("\r%d%%" %(float(self.i)/10))
			stdout.flush()
			
			if self.globalBest.fitness < goalFitness:
				break
				
	def particleSwarm2(self):
		v = 0
		self.globalBest = self.particles[randint(0, (len(self.particles)-1))]
		while self.i < max_iterations :
			for p in self.particles:
				fitness = Pso.circleFitness(p.params)
				if fitness < p.fitness:
					p.fitness = fitness
					p.best = p.params
				if fitness < self.globalBest.fitness:
					self.globalBest = p
				# v(t+1) = (w * v(t)) + (c1 * r1 * (p(t) - x(t)) + (c2 * r2 * (g(t) - x(t))) 
				v = (Pso.weight(self.i) * p.v) + c1 * random() * (p.best - p.params) \
						     		      	   + c2 * random() * (self.globalBest.params if self.i < 2 else self.findBestNeighbour(p).params - p.params)
				# x(t+1) = x(t) + v(t+1)
				p.params = p.params + v
			#saving the best Global Fitness per cicle for the plot
			self.fitnessVector += [(self.i+1, self.globalBest.fitness)]
			self.i  += 1
			stdout.write("\r%d%%" %(float(self.i)/10))
			stdout.flush()
			
			if self.globalBest.fitness < goalFitness:
				break

	def particleSwarm3(self):
		v = 0
		vmax = array([1]*2000) #max and min velocity
		vmin = -vmax
		repeatState = [0,0]
		self.globalBest = self.particles[randint(0, (len(self.particles)-1))]
		while self.i < max_iterations/2 :
			for p in self.particles:
				fitness = self.knapsackFitness(p.params)
				if fitness > p.fitness:
					p.fitness = fitness
					p.best = p.params
				if fitness > self.globalBest.fitness:
					self.globalBest = p
				# v(t+1) = v(t) + (c1 * r1 * (p(t) - x(t)) + (c2 * r2 * (g(t) - x(t))) 
				v = p.v + c1 * random() * (p.best - p.params) \
						+ c2 * random() * (self.globalBest.params - p.params)
				for j in range(len(v)):
					if v[j] > vmax[j]:
						v[j] = vmax[j]
					elif v[j] < vmin[j]:
						v[j] = vmin[j]
				# x(t+1) = x(t) + v(t+1)
				i = 0
				for par in p.params:
					if random() < 1/(1+exp(-v[i])):
						par = 1
					else:
						par = 0
					i += 1
			#saving the best Global Fitness per cicle for the plot
			self.fitnessVector += [(self.i+1, self.globalBest.fitness)]
			self.i  += 1
			stdout.write("\r%d%%" %(float(self.i)/10))
			stdout.flush()
			
			#if self.globalBest.fitness < goalFitness:
			#	break
			
	def particleSwarm4(self):
		v = 0
		vmax = array([1]*2000) #max and min velocity
		vmin = -vmax
		repeatState = [0,0]
		self.globalBest = self.particles[randint(0, (len(self.particles)-1))]
		while self.i < max_iterations/2 :
			for p in self.particles:
				fitness = self.knapsackFitness(p.params)
				if fitness > p.fitness:
					p.fitness = fitness
					p.best = p.params
				if fitness > self.globalBest.fitness:
					self.globalBest = p
				# v(t+1) = v(t) + (c1 * r1 * (p(t) - x(t)) + (c2 * r2 * (g(t) - x(t))) 
				v = p.v + c1 * random() * (p.best - p.params) \
						+ c2 * random() * (self.globalBest.params - p.params)
				for j in range(len(v)):
					if v[j] > vmax[j]:
						v[j] = vmax[j]
					elif v[j] < vmin[j]:
						v[j] = vmin[j]
				# x(t+1) = x(t) + v(t+1)
				i = 0
				for par in p.params:
					if random() < 1/(1+exp(-v[i])):
						par = 1
					else:
						par = 0
					i += 1
				
			#saving the best Global Fitness per cicle for the plot
			self.fitnessVector += [(self.i+1, self.globalBest.fitness)]
			self.i  += 1
			stdout.write("\r%d%%" %(float(self.i)/10))
			stdout.flush()
			
			#if self.globalBest.fitness < goalFitness:
			#	break
				
	def launchParticlesSwarm(self, problem):
		if problem == 1:
			return self.particleSwarm1()
		elif problem == 2:
			return self.particleSwarm2()
		elif problem == 3:
			return self.particleSwarm3()
		elif problem == 4:
			return self.particleSwarm4()
		
	def printParameters(self):
		print '\nPARAMETERS\n','-'*10
		print 'Population size       : ', self.population_size
		print 'Dimensions            : ', self.dimensions
		#print 'Goal Fitness	   : ', goalFitness
		print 'c1                    : ', self.c1
		print 'c2                    : ', self.c2
		print '-'*10
	 
	def printResults(self):
		print 'Global best fitness   : ', self.globalBest.fitness
		print 'Global best params    : ', self.globalBest.params
		print 'iterations            : ', self.i
		## Uncomment to print particles
		#for p in particles:
		#    print 'params: %s, fitness: %s, best: %s' % (p.params, p.fitness, p.best)
		
	'''
	def resetVariables():
		global i, globalBest, fitnessVector, particles
		i = 0
		particles = [] 
		globalBest = None
		fitnessVector = []
	'''
 
print '\nPARTICLE SWARM OPTIMIZATION\nby Marc Baro and Valerio Amodeo\n'
   
print 'Select the type of problem:\n1 - Circle Problem (fully connected)\n2 - Circle Problem (nearest neighbours)\n3 - Knapsack\n4 - Knapsack 0-1'
while True:
	try:
		problem = int (raw_input('Type a number: '))
		if problem > 4 or problem < 1:
			raise ValueError
		break
	except ValueError:
		print 'Invalid input - Not a number or wrong value'
if problem in (1,2):
	print 'Select the number of dimensions:'
	while True:
		try:
			dimensions = int (raw_input('Type a number: '))
			if dimensions > 2000 or dimensions < 1:
				raise ValueError
			break
		except ValueError:
			print 'Invalid input - Not a number or wrong value'
elif problem in (3,4):
	dimensions = 2000
print 'Select the number of particles:'
while True:
	try:
		population_size = int (raw_input('Type a number: '))
		if population_size > 1000 or population_size < 1:
			raise ValueError
		break
	except ValueError:
		print 'Invalid input - Not a number or wrong value'
print 'Select the run mode:\n1 - Normal\n2 - Statistics'
while True:
		try:
			mode = int (raw_input('Type a number: '))
			if mode > 2 or mode < 1:
				raise ValueError
			break
		except ValueError:
			print "Invalid input - Not a number or wrong value"
if mode == 2:
	print "How many runs?"
	while True:
		try:
			runs = int (raw_input('Type a number: '))
			if runs < 1 or runs > 1000 :
				raise ValueError
			break
		except ValueError:
			print "Invalid input - Not a number or wrong value"
c1 = uniform(1.0,2.0)
c2 = uniform(1.0,2.0)

if mode == 1:
	pso = Pso(population_size, dimensions, c1, c2)
	if problem in (1,2):
		pso.initParticles(True)
	elif problem in (3,4):
		pso.initParticles(False)	
		if problem == 3:
			pso.initPackages()
		elif problem == 4:
			pso.initPackagesVolume()
	pso.launchParticlesSwarm(problem)
	pso.printParameters()
	print 'RESULTS\n', '-'*10
	pso.printResults()
	print pso.valueWeight
	g = Gnuplot.Gnuplot()
	g.title('Best global fitness to iterations')
	g.xlabel('Iterations')
	g.ylabel('Global fitness')
	g.plot(pso.fitnessVector) 
elif mode == 2:
	k = 0
	pso = Pso(population_size, dimensions, c1, c2)
	pso.printParameters()
	print 'RESULTS\n', '-'*10
	while k < runs:
		pso = Pso(population_size, dimensions, c1, c2)
		if problem in (1,2):
			pso.initParticles(True)
		elif problem in (3,4):
			pso.initParticles(False)
			if problem == 3:
				pso.initPackages()
			elif problem == 4:
				pso.initPackagesVolume()
		pso.particlesSwarm()
		print 'RUN 			: ', k+1
		pso.printResults()
		k += 1